mdsquadfandomcom-20200213-history
Change position restraint force constant in MDP
By default, the force constant for position restraints in GROMACS is set in your topology. This makes it a pain to change. Martini has a good solution, which lets you set it in your MDP file. Here's how to implement that in any GROMACS system. Method A typical posre.itp file from pdb2gmx looks something like this: ; In this topology include file, you will find position restraint ; entries for all the heavy atoms in your original pdb file. ; This means that all the protons which were added by pdb2gmx are ; not restrained. [ position_restraints ] ; atom type fx fy fz 3 1 1000.000 1000.000 1000.000 7 1 1000.000 1000.000 1000.000 8 1 1000.000 1000.000 1000.000 9 1 1000.000 1000.000 1000.000 10 1 1000.000 1000.000 1000.000 Since GROMPP is derived from the C pre-processor, we can use some cleverness to replace this automatically. Typically we'll set the force constant in an mdp file, but if we don't then we want a default FC of 1000 units to apply (as it would in the above example): ; In this topology include file, you will find position restraint ; entries for all the heavy atoms in your original pdb file. ; This means that all the protons which were added by pdb2gmx are ; not restrained. #ifndef POSRES_FC #define POSRES_FC 1000.000000 #endif [ position_restraints ] ; atom type fx fy fz 3 1 POSRES_FC POSRES_FC POSRES_FC 7 1 POSRES_FC POSRES_FC POSRES_FC 8 1 POSRES_FC POSRES_FC POSRES_FC 9 1 POSRES_FC POSRES_FC POSRES_FC 10 1 POSRES_FC POSRES_FC POSRES_FC POSRES_FC is a variable. It can be defined in the MDP file; if it isn't, its set to 1000.00000 by default. Every time it occurs in the topology, it's replaced with the number when it's GROMPPed, and this gets you the correct restrained topology. A script to automatically process an ITP file is below. Once the topology is fixed, you just need to add a definition of POSRES_FC to your MDP. The MDP can set POSRES_FC with the define line. You'll still need to set POSRE, as otherwise the position restraint ITP won't be included and no restraints will be added - that is, setting POSRES_FC is not sufficient to turn on position restraints. For a force constant of 500, your MDP's define line should look something like this: define = -DPOSRES -DPOSRES_FC=500.0 Script The following is an AWK script to automatically process an ITP file. #! /usr/local/bin/awk -f # Josh Mitchell # July 2016 # Fixes a gmx pdb2gmx output posre_*.itp file so that the force constant can be defined in the .MDP with the line: # define = -DPOSRES -DPOSRES_FC=xxxxx # Initialise when we run into the position_restraints declaration /\[ *position_restraints *\]/ && fix_lines FALSE && init FALSE { init = TRUE # GROMACS' default force constant is 1000 default_fc = 1000 # We want to name the variable in the itp POSRES_FC fc_varname = "POSRES_FC" # Make it OK to fix individual lines fix_lines = TRUE # Stick a default force constant in if nothing is defined print "" print "; Modified by varposres.awk" print "" print "#ifndef", fc_varname printf "%s %s %4.6f \n", "#define", fc_varname, default_fc print "#endif" print "" # And then add in the declaration print $0 next } # If there's a second position_restraints declaration, crash /\[ *position_restraints *\]/ && init TRUE { print "Position_restraints declared twice" > "/dev/stderr" exit 1 } # If there's a new declaration, stop fixing stuff /\[ *A-Za-z0-9\_\-+ *\]/ { fix_lines = FALSE } # If $1 ~ /0-9+/ && $2 ~ /0-9+/ && $3 default_fc && $4 default_fc && $5 default_fc && fix_lines TRUE { printf "%5s %5s %10s %10s %10s \n", $1, $2, fc_varname, fc_varname, fc_varname next } # Print the line { print $0 } It will print the corrected ITP to STDOUT, so you'll need to redirect the output to a file. If you've saved it as varposres.awk in your $PATH, a command like this should work: mv posre.itp posre_bak.itp && varposres.awk posre_bak.itp > posre.itp Category:GROMACS Category:Guides